% 这个octave脚本基于达朗贝尔方程求解二次元电磁场
% 达朗贝尔方程将势写为波动形式，从而从势能导出场
% 由于没有正确设置边界，在边界附近的物理不可靠
% 其实我也不知道有没有计算错误...
%Gitee Repo

clc
clear

L=2;
dx = 0.01;
dt = 0.001;

% 波速
c = 1;

[x y] = meshgrid(0:dx:L);
n = size(x, 1);

% Amu=(\varphi,Ax,Ay,Az)是电势和磁势场
Amu0 = zeros(n,n,4);
Amu1 = zeros(n,n,4);
Amu2 = zeros(n,n,4);

sigma = (c*dt/dx).^2;

numSource = 1; %电流源个数，电流方向为+z方向
source(1).T = 0.5; %充电时间（的4倍）
source(1).x = 0.5*L; %位置
source(1).y = 0.5*L;

TICK=1000;
imgind=1;
for tick = 1:TICK
    t(tick) = tick*dt;

    % 达朗贝尔方程 1/c^2 \pdv[2]{A^\mu}{t} - \div A^\mu = \mu_0 j^\mu
    diff_i = Amu1(1:n-2,2:n-1,:)-2*Amu1(2:n-1,2:n-1,:) + Amu1(3:n,2:n-1,:);
    diff_j = Amu1(2:n-1, 1:n-2,:)-2*Amu1(2:n-1, 2:n-1,:) + Amu1(2:n-1, 3:n,:);
    Amu2(2:n-1,2:n-1,:) =2*Amu1(2:n-1,2:n-1,:) - Amu0(2:n-1,2:n-1,:) + sigma*(diff_i + diff_j);

    % 电流源
    for _n = 1:numSource
        _x = uint32((source(_n).x)/L*n);
        _y = uint32((source(_n).y)/L*n);

        if(_x < n && _y <n && _x > 0 && _y > 0) % 设置jz
            if tick*dt < source(_n).T/4
                Amu2(_y,_x,4) = sin(2*pi/source(_n).T*t(tick));
            else
                Amu2(_y,_x,4)=1; %恒定电流
            end
        end
    end

    %固定边界
    Amu2(1,:,:) = 0;
    Amu2(n,:,:) = 0;
    Amu2(:,1,:) = 0;
    Amu2(:,n,:) = 0;

	Amu0 = Amu1;
	Amu1 = Amu2;

    if mod(imgind,20) == 0
		clf
		hold on
		axis equal
        axis([0 L 0 L -1 1])
        xlabel('x')
        ylabel('y')
        caxis([-1 1])

        E = zeros(n,n,3);%Ex,Ey,Ez 是电场分量
        B = zeros(n,n,3);%Bx,By,Bz 是磁场（磁感应强度）分量

        %部分计算项没有被包括。

        %Ex = -\pdv{\varphi}{x} - \pdv{Ax}{t}
        E(2:n-1,2:n-1,1) = - (Amu1(2:n-1,3:n,1)-Amu1(2:n-1,1:n-2,1))/(2*dx) - (Amu2(2:n-1,2:n-1,2)-Amu0(2:n-1,2:n-1,2))/(2*dt);
        %Ey = -\pdv{\varphi}{y} - \pdv{Ay}{t}
        E(2:n-1,2:n-1,2) = - (Amu1(3:n,2:n-1)-Amu1(1:n-2,2:n-1,1))/(2*dx) - (Amu2(2:n-1,2:n-1,3)-Amu0(2:n-1,2:n-1,3))/(2*dt);
        % Ez = -\pdv{\varphi}{z} - \pdv{Az}{t}
        E(2:n-1,2:n-1,3) = - (Amu2(2:n-1,2:n-1,4)-Amu0(2:n-1,2:n-1,4))/(2*dt);

        %Bx = \pdv{Az}{y} - \pdv{Ay}{z}
        B(2:n-1,2:n-1,1) = (Amu1(3:n,2:n-1,4) -Amu1(1:n-2,2:n-1,4))/(2*dx);
        %By = \pdv{Ax}{z} - \pdv{Az}{x}
        B(2:n-1,2:n-1,2) = -(Amu1(2:n-1,3:n,4) -Amu1(2:n-1,1:n-2,4))/(2*dx);
        %Bz = \pdv{Ay}{x} - \pdv{Ax}{y}
        B(2:n-1,2:n-1,3) = (Amu1(2:n-1,3:n,3) -Amu1(2:n-1,1:n-2,3))/(2*dx) - (Amu1(3:n,2:n-1,2) -Amu1(1:n-2,2:n-1,2))/(2*dx);

        view(30,60);
        range = 1:10:n;
        zeromat = zeros(size(x(range,range)));

        quiver3(x(range,range),y(range,range),zeromat,E(range,range,1),E(range,range,2),E(range,range,3));
        quiver3(x(range,range),y(range,range),zeromat,B(range,range,1),B(range,range,2),B(range,range,3));
%        quiver3(x(range,range),y(range,range),zeromat,Amu1(range,range,2),Amu1(range,range,3),Amu1(range,range,4));

		drawnow
		pause(0.01)
    end
    imgind++;
end


